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TREATMENT OF INCOMPATIBLE INITIAL AND 
BOUNDARY DATA FOR PARABOLIC EQUATIONS IN 
HIGHER DIMENSION 

QINGSHAN CHEN, ZHEN QIN, AND ROGER TEMAM 



Abstract. A new method is proposed to improve the numeri- 
cal simulation of time dependent problems when the initial and 
boundary data are not compatible. Unlike earlier methods limited 
to space dimension one, this method can be used for any space 
dimension. When both methods are applicable (in space dimen- 
sion one), the improvements in precision are comparable, but the 
method proposed here is not restricted by dimension. 

1. Introduction 

When performing large scale numerical simulations for evolutionary 
problems, we use most often initial and boundary conditions provided 
by approximations, by other simulations, or by experimental measure- 
ments. These data may not satisfy certain compatibility conditions 
verified by the solutions; thus various modifications deemed non essen- 
tial are made on the data to overcome these difficulties. Such issues 
are extensively addressed in the literature; see for instance in geophys- 
ical fluid mechanics [3] or [27] which contains many allusions to this 
difficulty; in classical fluid mechanics, see e.g. [71 [121 HSl HSl US]; see 
also [1] in chemistry and [28] in a general mathematical context. 

We want to address here a less known difficulty of "mathematical" 
nature which, the specialists believe, will become very important as we 
move to high resolution methods thanks to the increase of computing 
power and memory capacity of the computers. A very simple exam- 
ple of such a difficulty appears when solving in space dimension one 
on (0,1), the heat equation ut — u^x = with boundary conditions 
u{0,t) = u{l,t) = and initial condition u{x, 0) = 1. The solution ex- 
ists and is unique (for t > 0) and the analytic expressions of u are pro- 
vided in the literature (see e.g. [1]). This problem is simple enough that 
it can be solved satisfactorily by numerical methods, but the solution 
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does display singularities in the corner x = 0, t = and x = l,t = 0. 
For this problem and for general parabolic equations, it is known from 
semi-group theory [HI [22] or by using the analyticity in time of the 
solutions (see [11]) that certain norms of du/dt grow as a power of 1/t 
when t — 7- 0. It is believed that such singularities will affect large scale 
computations as our demand for better results increases. In fact it 
has been observed by some authors that, when using spectral methods 
for the space discretization, the spectral accuracy is lost if nothing is 
made to address this singularity and a series of works resulted from this 
observation; see e.g. [21 El El 13 [ID] , see also [S] for a nonlinear equation. 

The mathematical difficulty studied in detail in e.g. [IHl EHl [2D1 1231 
1211 125] is the following; even if the initial and boundary data of an evo- 
lution problem are given C°°, the solution may not be C°° near t = 0. 
In fact, k compatibility conditions between the data are needed for the 
solution to be near t = and hence an infinite number of compati- 
bility conditions are needed for the solution to be C°° . Furthermore the 
initial and boundary conditions that are compatible form a relatively 
small set (in an informal sense), so that most numerical simulations 
are done with data which are not compatible, generating a loss of ac- 
curacy near t = if nothing is done. In the works mentioned above, 
methods have been proposed to address the first or the first two incom- 
patibilities. It is believed that dealing with one or two incompatibilities 
substantially improve the quality of the simulation and, in any case, 
dealing with more incompatibility conditions may become impractical. 
However a strict limitation of these works is that the proposed methods 
only apply to space dimension one and, to the best of our knowledge, 
there is (there was) no method available in dimension two or larger to 
address this difficulty. 

As we said, past works, on the computational side, have been de- 
voted to space dimension one. This problem has been addressed in a 
series of articles by Flyer, Boyd, Fornberg and Swarztrauber [21 UHl El [9] 
who proposed a number of remedies in space dimension one for linear 
equations. Nonlinear equations in space dimension one were considered 
in [g. 

In these articles, the authors introduce a correction term in the linear 
and nonlinear cases by setting 




u = V + S, 
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where 5* absorbs the incompatibihties between the initial and bound- 
ary data up to a certain order. Now, free of incompatibihties of lower 
orders (the most severe ones) , v is computed by an appropriate numer- 
ical procedure, such as finite differences, the Galerkin finite element 
method, spectral or pseudo-spectral methods. As a final step, the 
original solution u is recovered through (11.11) . This remedy procedure 
effectively reduces the errors at the spatio-temporal corners during the 
short initial transient period. 

For dimensions higher than one, the construction of 5*, to correct for 
singularities generated at t = by incompatible data, remains an open 
problem. A method to overcome this difficulty is proposed, analyzed 
and tested in this article. 

In this article, we intend to study, from both a theoretical and numer- 
ical point of view, the incompatibility issue for the multi-dimensional 
time-dependent linear parabolic equation: 



We believe that our method applies to more general parabolic equa- 
tion, but we restrict ourselves to equation (II. 2p in this article devoted 
to feasibility. 

The method that we propose is based on the concept of penalty. 
We replace in (II. 2p the boundary value u\Qn = 9 by u\do. = k'^ ■ This 
boundary value fc^ which depends on a parameter e > is such that 
k^\t=o = Mo Ian (see equations (12.11) . (12. 2p below), so that the first incom- 
patibility has disappeared. Now /c^, through a penalty procedure with 
parameter e, is forced to rapidly vary from Mo|an at t = to the desired 
value, namely g. This is achieved through equation (12.21) : initially kl is 
large, but it becomes rapidly of order 1, and then, by the first equation 
(12.21) . A;^ — (yf is of order e. It is easy of course to integrate equation 
(12.21) although an explicit solution is not available in the general case 
where g depends on time. The concept of penalty has been introduced 
in the mathematical literature by R. Courant [6J; it has been adapted 
to evolution problems by J. L. Lions in [21j, a reference which contains 
many evolution equations similar to (12. 2p (Chapter 3, Sections 5 to 8); 



(1.2) 




4 QINGSHAN CHEN, ZHEN QIN, AND ROGER TEMAM 

it is widely used in optimization 0; see also [26] (Chapter 1, Section 
6). In this work, we firstly present our approach in details and study 
it theoretically to prove the strong convergence of the method. Then 
we implement it numerically on a number of examples. Because the 
penalty method does not depend on the properties of fl, we believe 
that this method can be applied to many systems with many different 
domains Q. The question that remains is the choice of e > small. In 
optimization theory, the choice of e is usually made by trial and error 
and is not a major issue. It does not follow the "intuitive" idea that 
the error becomes smaller as e becomes smaller because of many other 
contingent errors such as round-off and descretization errors. In general 
the error becomes "optimal" for some value of e and the method gives 
less good results for smaller or larger values of e. In our case (see Fig. 
inland [7]), at the initial steps, the error decreases sharply as e increases 
and remains close to 0, then it becomes stable flat. At the final steps, 
the error increases almost linearly as e increases. With e at about 0.1, 
the initial error is minimized while the error at the final step is well 
controlled. In a short time period e = 0.5 gives us smaller errors and 
again after a short time period e = 0.1 gives us a smaller errors. In 
general the choice of e really depends on our goals of the computation. 

This article is organized as follows. In Section 2 we present the 
method and establish various approximation results. Then, in Section 
3 we present numerical results showing the efficiency of the method 
and comparing it to earlier methods. In Section 4 we present some 
conclusions and perspective of future developments. 



2. PENALTY METHOD 

2.1. Perturbed problem (and the statement of the main re- 
sult). We consider the system (11. 2p . where z/ > 0. If Uo\qq ^ fl'(O), 
then we face an incompatibility problem, in which case we consider a 
new system instead, namely, for e > fixed, 

u'\t=o = Mo, 



search on Google with the words "optimization, penalty" produced 3,350,000 
entries. 
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^22) kt + \k^~g)=0, teR^, 

k'{0) = uolan- 

In this article, | ■ | is the L'^{Q) norm, and || ■ ||= |V ■ | is the Hq{Q) 
norm; for other norms, we will use the subscript notation. 

The system fl2.ip - fl2.2l) is actually decoupled and fl2.2p is just an 
Ordinary Differential Equation with x G dQ as a parameter. As we 

see below, if we are given g, g' = & L'^{0,T; H^{T)), then we 

have the existence and uniqueness of k'^ in L^(0, T; ifa (F)) and fur- 
thermore, by the effect of the penalty term, {k"^ — g)/e,k^ converges to 
g in suitable spaces as e — >■ 0. Equations (12.11) is a heat equation with 
non-homogeneous boundary conditions, and we have the existence and 
uniqueness of a solution if the data are sufficiently regular. Then we 
have the following theorem. 

Theorem 2.1. Assume that we are given g G L°°{0,T; H^{T)) (F = 
dn), with gt e L'^{0,T;H^T)), and Uq G H\n). Then (EJ) has a 
unique solution u G L^(0, T; i^Q fl C([0, T]; and for each 

e > 0, /fO)-/fOI) has a unique solution G L\0,T; H\n))n C([0,T]; 
L^{n)), k' G L2(0,T;if5(r)). Furthermore, ase ^0, 



^ u in L^{0,T;H ^(i^)) strongly, and in 
^^'^^ C{[to,T];H-^{n)) strongly, Vto > 0. 



Remark 2.1. We do not prove a strong convergence of to u on all 
of [0,T] in the L°° sense, and we do not expect such a convergence 
to occur since u has a singularity at t = 0. Alternatively one could 
capture the singularity of u near t = by using the methods of singular 
perturbation theory as in, e.g. Jung-Temam [T7], which we do briefly 
in Section 12. 3^ and will also be studied elsewhere. 

Before we prove Theorem 2.1, we will first prove the following lemma. 

Lemma 2.1. If g e L^{0,T; H^T)) and g' G ^^(0, T; ifl (E)), 
then there exists a unique k^ m L^{0,T;H2{T)) satisfying / fOI) . and 
as e ^ 0, k^ g in L^(0,T; H^iT)) strongly. Furthermore, as e 0, 

/ k%s)ds^ / g{s)ds m L'^{0,T;H^T)) strongly. 

Jo Js 
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Proof. We explicitly solve the ODE system f l2.2p . and we obtain the 
solution F e L2(0,T; H^T)): 

(2.4) F(t) = e-^F(O) + / -g{s)e'^ds. 



e 



Then we rewrite ( I2.2p i in the form 



(2.5) (F-^), + i(F-^7) = -^7,. 



Taking the scalar product of (12. 5 p with k"^ — g in (F), we obtain 



< IqA 1 \k^ — q\ 1 

- 2' 'H^(r) 2£' ^'H3(r) 



Hence 



(2.6) ^\k'-g\\ +-\k'-g\\ < e\gt\'^ i . 
Using the Gronwall inequality we obtain 

(2.7) |A;"-(?ri it) < e-"A¥ - g? . (0)+£|(7tP i • 
We integrate (EZD over (0,T), and we obtain (mq G ^T^^) @) 

(2i 



T 







|A;'-^?l^(^^^^^<^(l-e--)Mo|^-^7(o) 



2 



2 



1 - , 

H^r) ' 'L2(0,T;H^(r))' 



and hence 



which implies 

(2.9) k' ^ g strongly in L'^{0,T; H^T)) as e ^ 0. 



We do not address the question of minimal regularity of uq, that is e.g. uq G 
i^(f2), which is not in the scope of this article. Indeed the problem of incompatible 
data occurs already with very smooth data. 
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Integrating fl2.5p from to t, we obtain 




which yields 




Jo 



The proof of Lemma 12.11 is complete. 



□ 



Remark 2.2. We could prove a stronger result namely, k^ g, / k^{s)ds 



/ g{s)ds strongly in L'?(0, T; iJ2 (r)), for all 1 < g < oo. But in 



this article, g = 2 is enough for our needs. And for q = oo, from (12. 7p . 
we obtain 

(2.12) k'-g = 0{V^) in L°°(to,T;i75(r)) for Vto > 0, 

and also from (12.101) . because k^—g and g are bounded in L~(0, T; H2{r)), 
we obtain 

(2.13) [ k%s)ds^ [ g{s)ds strongly in L°°{0,T;H^{r)), 



the norm of the difference being of order e. 

2.2. Convergence results for u^. Since Q is smooth, there exists a 
lifting operator L, linear continuous from H2(r) to H^{Q). We con- 
sider such an operator and set K"^ = L{k'^),G = L{g), and thus have 



by assumption G G L°°(0, T; ffi(fi)), Gt G L'^{0,T; H\n)). So we 
immediately infer from dMD, f l^TTT]) . fl^TT^ and fin^ that, as £ ^ 
(2.14) 

^G strongly %n L\0, T; H\n)) n L°°(to, T; H\n)) Vto > 0, 



We now prove Theorem 12. 1[ 

Proof of Theorem 12.11 Set t>^ = — K'^; then the system (12.11) 





(2.15) 




yields 
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(2.16) 



vl - uAv' = f + uAK', 
f^|t=o = Mo - ^^(0) = Mo - Luolan, 
v'ldu = 0. 



Integrating (l2.16P i from to t, we obtain 
(2.17) 



v%t)-uA [ v%s)ds= [ f{s)ds - K%t) + uA [ K%s)ds + uo. 
Jo Jo Jo 

We set V = [ v'{s)ds (withF"(0) = 0), and F{t) = [ f{s)ds+Uo; 
Jo Jo 



so solves the following system 
(2.18 



Vf - vAV = F - + vA [ K%s)ds, 

Jo 



V'\t=o = 0, 

We take the scalar product of (12. 181) i with in L^(f2) and find, 
(2.19) 

~\V'\^ + iy\\V' f = (F, V) - (K^ V) + u{A [ K'{s)ds, V). 
2 dt Jq 

We can bound the terms in the right-hand-side of (12.191) as follows: 

(2.20) {F,V') < \F\\V'\ < Ci\F\ \\ V \\< c[\F\^ + - || V f , 

6 

(2.21) - (K',V') < \K'\\V'\ < ci\K'\ II V \\< c' liT^I^ + - || V f , 

6 

z/(A [ K'{s)ds,V') = [ K%s)ds,VV') 

Jo Jo 



(2.22) 



< z/ II / K'{s)ds nil V 

'0 



■3 II / K%s)ds f . 
Jo 



Here and below, the c,c',Ci,c[ are various constants independent of 
e, which may be different at different places. 
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Combining fICTD . (^Mj, (12:2T1) and ( Km gives 

(2.23) -\V'\^ + u\\V' \\^<c[\F\^ + c'^\K'\^ + 4\\ / K%s)ds f 
at Jq 

Integrating fl2.23p over (0,t), we obtain 



t 

e\2 



(2.24) 



\V'{t)\^ + u [ \\V'fds<c[[ \F\^ds + d^f \K' 
Jo Jo Jo 

+ 4 [ W [ K'{T)dT IP ds. 

Jo Jo 

We also integrate f l2.23p over (0,T), and obtain 



ds 



T 

'^^ds 



\V'{T)\^ + u [ \\V'fds<c[[ \F\''ds + d^f \K' 

Jo Jo Jo 

+ 4 /" W [ K%T)dT f ds. 

Jo Jo 



(2.25) 



It follows from fl^l^ and KWi that and j K'{s)ds are bounded 

in L'^{0,T;H\Q)), and thus (ICTD . (12:251) yields: 
(2.26) 

remains bounded in L°^{0,T; L'^{n))nL^{0,T; H]{n)) as e ^ 0. 

Thus, there exists a subsequence V' and 1/ G L°°(0, T; L2(fi))nL2(0, T; 
ifQ(f2)) such that, as e' — )■ 0, 

(2.27) V"' ^ V weakly in L'^{0,T; H^{n)), 

and weak — star in L°°{0,T] L'^{Q)). 

Using fl2Ail) . fl2:T5|) and fl2:27D . we can pass to the limit in fl236D 
with the sequence e' — 0. We proceed as follows. 

For all a G H^{^), and </) in Ci(0,T) with 0(T) = 0, we multiply 
(l2.18P i by acj) and integrate over Q x (0, T); we obtain 
(2.28) 

- /" {V'',a)(f)'{t)dt + iy [ {VV'\Va)(t){t)dt= I {F,a)(pit)dt 
Jo Jo Jo 

{K'\a)(t){t)dt-v I {V f K''{s)ds,Va)(p{t)dt. 
Jo Jo 
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Passing to the limit with fl214ll . (EJS]), (ESZD, we find 
(2.29) 

{V,a)(p'{t)dt + u I {W,Va)(p{t)dt= [ {F,a)(P{t)dt 



{G,a)(p{t)dt-u / (V /" G{s)ds,Va)^{t)dt. 
Jo Jo 

Taking G V{0,T), we see that V satisfies 
(2.30) 

{Vt,a)-u{AV,a) = {F-G + iyA [ G{s)ds,a), V a e ifo^(fi). 

Jo 

Now we want to show that V^(0) = 0. 

We classically integrate fl2.30p times 0(t) over (0, T) and we obtain: 
(2.31) 

- / {V,a)(p'{t)dt + u [ (VV, Va)0(t)rft = / (z/A [ G{s)ds, a)(p{t)dt 
Jo Jo Jo Jo 

+ [ (F-G,a)0(t)c/t + (V(O),a)0(O). 
Jo 

By comparing with (I2.29p . we find that 
(2.32) (V(O),a)0(O) = O, 

for every a G Hq{Q) and every G C"'^([0,T]) with 0(T) = 0. This 
implies ^(0) = as desired. Finally V satisfies 



(2.33) 



Vt-uAV = F -G + uA [ G{s)ds, 

Jo 



V\t=o = 0, 
Vlan = 0. 



Remark 2.3. Furthermore, we could prove that the whole sequence 
V weakly in L'^{0,T; H^{n)), and weak star in L°°(0,T; L'^{n)). 
Indeed, if not, arguing by contradiction, we could find a subsequence 
— )■ 0, such that 



V'^ in L\0,T;H^{Q)) weakly, 

L°°(0,^;L2(^])) weak -star. 



Repeating the argument above leading to (12.271) . we could extract from 
Ei a subsequence e[ and find V such that, as e[ — >• 0, 
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V< in L\0,T;H^{Q)) weakly, 

L°°(0,T;L2(n)) weak -star, 



(2-35) 



where V is the solution of (I2.33p . But the solution of fl2.33p is unique; 
hence V = V, and then (12.351) contradicts (12 .34^ . 

Before we finish the proof of the theorem, we now prove the following 
Lemma. 

Lemma 2.2 Under the assumptions of Theorem \2.1\ with V, being 
the solutions of h2. 33\) and l\2.18\) . we have, as e 0, 

(2.36) V strongly in L\0, T; H^{n)) n C([0, T]; L\n)). 

Proof. We subtract ([2718]) i from f l233|) i. and obtain 
(2.37) 

Vt-Vt'-u{AV-AV') = K'-G + u{A [ G{s)ds-A [ K'{s)ds). 

Jo Jo 

We then take the scalar product of (I2.37P with V — V^ in L'^{fl), and 
integrate in time from to t, and we obtain: 
(2.38) 



.\iv-v^m\' + u 



■[ \\V -V f ds= [ {K' -G,V -V')ds 
Jo Jo 
t pt pt 

(z/A / G{s)ds-vA I K'{s)ds,V -V')ds. 



10 Jo Jo 
Now we set 

"t 



(2.39) Xe{t) = \\{V-V%t)\'+''-j^ \\V-V^rds, 

and estimate the right-hand-side of (12.381) as follows: 
(2.40) 

f {K' -G,V -V')ds< [ \K' - G\\V - V'\ds 
Jo Jo 



<ci [ \K' -G\\\V -V II ds 
Jo 

<ci{! \K' -G\^dsy^{ [ \\V-V'f 
Jo Jo 

<c I \K'-G\^ds + ^f \\V-V'fds 



ds)' 
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(2.41)^ 

z/(A / G{T)dT-/\ / K'{T)dT,V -V')ds 



<^ II / {Gi.^)- K'{T))dT WW V -V \\ds 

Jo Jo 

<c' [ W [ iG-K')iT)drW^ ds+'^ [ WV-V'W^ds. 

Jo Jo 4 Jq 

Combining ( 12.40p and f l2.4ip . we see that 

(2.42) xe{t) < c' r W [\g - K'){T)dr f ds + c f \K' - G^ds 

Jo Jo Jo 

The right-hand side of (12 .42^ converges to as e converges to 0, 
because of (I2.14p and (12.151) . and so does Xeij)- For t = T, we find 

(2.43) V strongly in L^{0,T; H^{n)) as e ^ 0, 
and taking the supreme of (12.421) with respect to t, we see that 

(2.44) V strongly in L°°{0,T; L\n)) as e ^ 0. 

The lemma is proved. □ 

Now we apply Lemma [2.21 and obtain as e — )■ 0, 
(2.45) 

AV -^AV strongly in L\0,T; H-\Q)) nC{[0,T]; H-\Q)), 

and from (I2.15p . we obtain as £ — )■ 0, 
(2.46) 

A [ K%s)ds^A [ G{s)ds strongly in L'^{0,T; H-\n)), 
Jo Jo 

so after comparing (12. 181) i with (l2.33P i . we conclude that as e — )■ 0, 

V/ ^Vt strongly in L'^{0,T; H-\n)) and 
^^■^'^^ C{[to,T];H-') for Vto > 0. 



Now we define v = Vt, and take the derivative on ( 12.33^ 1 . we obtain 
that V solves the following system, 

Vt - vAv = f -Gt + uAG, 
(2.48) { v\t=o = uo-G{0), 

v\dn = 0. 
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So (12:471) yields as e ^ 0, 

^ V strongly in L'^{0,T; H~^{il)) and 
^^■^^^ C{[to,T];H-') far Vto > 0. 



The final stage of the proof of Theorem 12.11 consists in reinterpreting 
the results above that is (12.141) . fl2.15p and fl2.49p in terms of the con- 
vergence of = f ^ + towards u = v + G; we obtain precisely (12. 3p . 
Theorem 12.1! is proven. □ 

Remark 2.4. Similarly as Remark 12.21 we see that we also have ^ u 
strongly in Li{0,T; for all 1 < g < oo. 

2.3. Boundary layer analysis for /c^. In the previous section, Lemma 
12.1! stated that under our assumptions, strongly converges to g in 
L^(0,T; if2(r)), as £ — 0. Here in order to better compare k'^ and g, 
we are going to study the boundary layer for the system (12. 2p . 

Along the asymptotic analysis, we define the outer expansion k^ ~ 

oo 

^^e-' A;-' . By formal identification at each power of e, we obtain 

j=0 

0(6-') : k' = g, 

(2.50) 

0(£^) : kl + = 0, Vj > 0. 
By explicit calculations, we find: 

(2.51) y = {-iyg^^\ Vj > 0. 

It is clear that the functions k^ of the outer expansion do not gener- 
ally satisfy the initial condition in (12. 2p in the case of interest here where 
5^(0) UqIqq. To account for this discrepancy, we classically introduce 
the inner expansion k^ ~ Yl'jLo^^^'' y where 9^ = 9^{t), (t = t/e). Then 
we find 

oo j oo 

By formal identification at each power of e, we obtain the following 
equations: 

(2.52) +^^(t) = o, /or j>0. 

at 
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The initial conditions that we choose are: 

e\0) = uo\9n-g{0), 
' 0^(0) = -y{0) = {-iy+^g^^\0), for j > 1. 

By exphcit calculations, we obtain: 

(2.54) = e-ie^{0), for j > 0. 
To obtain the asymptotic error estimate, we set 

(2.55) Wgn = k ksn ^en; 

where 

n n 

ken = ^ ^ k^ , 9f,ri = ^ ^ 9^ . 
j=0 j=0 

Now we can conclude as follows. 

Theorem 2.2 // G L'^{0,T; H^T)) for n > 0, and w,n is 

defined in ^2. 55]) . then as e — )■ 

Wsn = 0{E^^^') m L\Q,T-hHt)), 

(2.56) , , 
^i;,„ = 0(£"+2) L°°(0,T;i72(r)). 

Proof. We firstly notice that w^n vanishes at t = 0. We insert then 
fl235|) into (D, and we find: 

Wen\t=Q = 0. 



(2.57) 



We take the H^iT) scalar product of f l2.57l) i with Wen and integrate 
over [0,t]; we obtain 

-\Wenit)\\ +t\w,Js)\\ ds= l\{-eY+^g^''+^\wen) ^, 4S 
1 /■* p2(n+l) rt 

(2.58) 

£k.n(t)|'i + /" \Wen{s)\\ < e^C^+l) / | (s) rfs. 

' ^'//^(r) ' ^'//5(r) - Jo ^'H5(r) 

If we set t = T in f l238l) . we obtain Wen = 0(£"+i) in L2(0, T; //^(r)), 
and if we take the supremum of (12.581) over [0,T], we obtain Wsn = 
0(e"+^) in L°°(0,T;i73(r)). 

Theorem 12.31 has been proved. □ 



TREATMENT OF INCOMPATIBILITIES 15 

Remark 2.5. If we additionally assume that gtt e L^{0,T;H2) in The- 
orem [2TT], from fl2.56p . setting n = 1, we find 

Wei = k'-g + tgt - {u^\dn - 9me~"' + £^t(0)e-*/^ = ©(^^/S), 
in L°°(0,r;iy^(r)). Then for any to > 0, 

{k' - g){to) = -egtito) + 0{e^/^) + e.s.t., 
where e.s.t. means exponentially small term (for all ilf^-norms). 

Hence Kk"^ ~g){to)\^i = 0{e), for Vto > fixed. We now do similarly 
as fl2.5l) - fl2.7p for kf — gt and integrate from to to t, 

(2.59) |M-c/tpi it) < e-'-^M - gA^ ^ {to) + e\gu? i , 
which yields 

(2.60) kl-gt = 0{y/e) in L°°{to,T; H^T)), Vto > 0. 

3. Numerical Results For the penalty method 

3.1. Approximations of k"^. In order to test the efficiency of the 
proposed penalty method, we will provide in this section and in the 
next one (Sec. 13.31) some numerical results for system (12. 2 p with Q = 
(0, 1) X (0, 1), and < t < 1; we set g{t) = sin{t) for all {x,y) e dQ, 
UQ{x,y) = sini^x + ^)sm(^y + ^), in which case, we face the dis- 
crepancies all along the lines x = and y = (but no discrepancy 
along the parts x = 1 oi y = 1 oi the boundary). 

We start by testing the quality of the approximation of fc^ inferred by 

n 

the boundary layer analysis of Section [231 that is k^ ~ e^{y + 9^)^ 

j=0 

for suitable e's and n's. 

Because (12. 2p is a 2D system which makes the graphing impossible 
along the time axis, we restrict ourselves to follow the time evolution of 
the exact and approximate function at one point of the boundary; for 
simplicity we choose the corner (x, y) = (0, 0). We then plot in Figure 

n 

m g{t) (solid line) and F (dash-dot line), k' ~ J]£^(fc^' + 6^) with 

j=0 

n = or 1, and e = 0.1 or 0.01. For n = 1, the proposed new scheme 
gives a good approximation of fc^. 
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n=0 and epsilon=0.1 
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(a) Time 

n=0 and epsilon=0.01 
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Figure 1. (a) Boundary Layer Element with e = 
0.1, /c*^ ~ A;'^ + (P ^ (b) Boundary Layer Element with 
e = 0.1, k° + 9° + ek^ + eO^, (c) Boundary Layer 
Element with e = 0.01, k^ ~ k^ + 6^, (d) Boundary Layer 
Element with e = 0.01, fc^ ~ + ^° + ek^ + e9^. 



Fig. |2] gives the L^— and L°°— errors which stand respectively for 
the L\0,T;H^r)) and L°°{0,T] H^T)) norms of the difference be- 

n 

tween the real solution k'^ and the approximations ''^^e^^y + 9^), as 

j=0 

the number of time steps T, e and n vary. It is clear that the smaller 
e is, the smaller both errors are. 



3.2. A two-dimensional system in a square Q. To verify the effec- 
tiveness of the penalty method, we use the finite elements method for 
the spatial approximation of u. The Penalty Method is mainly aimed 
for multi-dimensional time-dependent PDEs, so we consider the 2D 
system, as in (11.21) : 
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n=0,T=10 




-4 -3 -2 -1 

(c) Logarithm of Epsilon (Base 10) 



n=0,T=100 
S- 0.5 1 • 




(b) Logarithm of Epsilon (Base 10) 



n=1,T=100 




-4 -3 -2 



(d)Logarithm of Epsilon (Base 10) 



(3.i; 



Figure 2. ^^(plot o)- and L°°(plot *)- error between 
the boundary layer schemes and the real solution for ekt+ 
k = g. 



— - V[U^^ + Uyy) = /, 



u\dn = 9, 
u\t=o = Uq. 

where 0<a;<l,0<?/<l,0<t<l. 

Stt Stt Stt Stt 
We set z/ = 0.2, / = 0, g = 0, Uq = sin{—x + —)sin{—y + — ), 

so that g{0) ^ uo\gn on the lines a; = and y = 0. So we face the in- 
compatibility problem, namely the boundary conditions and the initial 
condition do not match at these corners of the time and spatial axes. 
For the test we set e = 0.1 in the penalty approximation fl2.ip -f l2^ 
of (13. ip . In more general problems, we might have discontinuities at 
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solution at t=0 solution at t=0.5 solution at t=1 




Figure 3. The exact solution of the system in the 
square Q without applying the penalty method, at times 
0,0.5 and 1 (Figures [3(a), El^b), El^c)). 



section of the solution at y=0.6,t=0 section of tfie solution at y=0.6,t=0.02 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



X X 



section of tfie soiution at y=0.6,t=0.04 section of the solution at y=0.6,t=0.08 




"0 0.2 0.4 0.6 0.8 1 "0 0.2 0.4 0.6 0.8 1 



Figure 4. The sections of the exact solution at y = 0.6 
when t is close to 0. 

the space corners x = 0orl,?/ = 0orl. But since the function uq is 
smooth at the corners, these singularities do not occur here, at least at 
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Figure 5. The comparative errors of the 2D system in 
the L°° norm for e = 0.1. (a) Maximum comparative 
error for a short time period (in real value), (b) Maximum 
comparative error when applying the Penalty Method 
[times 10~^). 

the low orders. 

We first plot the solution of system fl3.ip without applying the penalty 
method. The solution is plotted in Fig. El (a) is the graph of the 
approximate solution at t = 0, (b) is the graph of the approximate 
solution at if: = 0.5, and (c) is the graph of the approximate solution 
at t = 1. The graph displays a sharp gradient around the corner of 
the time-space axis during an initial short period due to the incom- 
patibility between the initial and boundary conditions there. In order 
to see the sharp gradient clearly and the changes of the gradient as 
time evolves, we plot the sections (x G (0, = 0.6) of the solution 
at times close to 0; see Figure ID It is clear that, at t = 0, we observe 
the sharpest gradient at the time-space corner and as time evolves, the 
gradient becomes smoother and smoother at that corner, until t = 0.08, 
when it is essentially fiat. 

Next, to study the accuracy of the numerical method, we must mea- 
sure the errors for the approximate solutions. Hence we compute the 
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Penalty Error at Initial Steps Penalty Error at Final Steps 




0.5 1 0.5 1 

Epsilon Epsilon 



Figure 6. The comparative errors for the 2D system 
in the norm with variations of e with mesh Ax = 
^, Ay = j^, At = (a) The error at initial step, (b) 
The error at final step. Note the factors 10~^, 10~^ in 
(a), (b) 




Figure 7. The maximum comparative errors for the 2D 
system in square domain. 



comparative errors which are the differences between two numerical so- 
lutions for the problem, one with the stated mesh sizes, and the other 
one with a finer mesh. Then at each time step, we obtain the maxi- 
mum error between the two meshes above; it is understood to be L°° 
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Comparative Maximum Error at Eacti Time Step 




Figure 8. The comparative errors for the 2D system 
in sense at e = 0.1. The upper hne is with mesh 



Ax 



mesh Ax 



24' 



Ay 



48 



24' 



Ay 



At 

48' 



1000' 

At- 1 



the lower hne is with 



4000 



comparative errors, or maximum comparative errors. In what foUows, 
aU the error terms are to be understood in this sense. 



We plot the maximum comparative errors of the 2D system on Fig. [51 
Graph (b) is the plot of the maximum comparative errors along the 
whole time period if we apply the penalty method. Because the dis- 
crepancy happens at the time-space corner, we zoom into the left corner 
of graph (b) and compare it with the error when we do no apply the 
penalty method. In graph (a), the line with stars is the maximum 
comparative errors with the penalty method applied, and the line with 
circles is the maximum comparative errors without the penalty method. 
We observe that the magnitude of the errors at the time-space corner 
is reduced by around one order by the penalty method. 

Because we use finite differences method, so for the same e, if we 
have a finer mesh, the maximum comparative error should be smaller. 
In Fig. [H]we plot the error of the 2D system with e = 0.1, the lower 
curve with a finer mesh, the upper curve with a coarser mesh. The 
magnitude of the errors are reduced by around 40%. So for a fixed e, 
the finer the mesh is, the smaller the error is. We are also interested in 
the decay of the maximum errors. The most interesting and informa- 
tive comparison can be made between the decay rates of the maximum 
errors at the initial and final time steps. In Fig. [9] (a), the maximum 
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Maximum Error at the Initial Time Steps Maximum Error at tine Final Time Steps 



— * — Without Penalty Method 
4K ^< * * — e — With Penalty Method 




logarithm of the number of the logarithm of the number of the 

(a) segment in X (base 10) (b) segment in x (base 10) 



Figure 9. Decay of the maximum errors. When we 
apply the penalty method here e = 0.1. (a) at the initial 
steps, (b) at the final steps 



errors at the initial time step are plotted against the grid resolution 
in the loglog scale. Without the penalty method, the maximum er- 
rors do not decrease as the grid refines, which demonstrates that the 
singularity in the solution during the initial period is serious. With 
the penalty method {e = 0.1), the maximum errors decay at roughly 
the second order. Fig. |9] (b) shows that, with and without applying 
penalty method, the maximum errors at the final steps (t=l) decay at 
approximately the second order. 

Next, we fix the meshes at e.g. Ax = ^, Ay = ^, At = and let 
e vary. In Fig. |6l we plot the maximum comparative errors of system 
(13. ip . At the initial steps, the error decreases sharply as e increases and 
remains close to 0, then it becomes stable fiat. At the final steps, the 
error increases almost linearly as e increases. With e at about 0.1, the 
initial error is minimized while the error at final step is well controlled. 
But as Fig. [7] shows, in a short time period e = 0.5 gives us smaller 
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solution at t=0 solution at t=0.5 solution at t=1 




(a) (b) (c) 



Figure 10. The solution of the system in disk f2 with- 
out applying the penalty method. 



errors and again after a short time period e = 0.1 gives us a smaller 
error. In optimization theory, the choice of e is usually made by trial 
and error and is not a major issue. It does not follow the "intuitive" 
idea that the error becomes smaller as e becomes smaller because of 
many other contingent errors such as round-off and descretization er- 
rors. In general the choice of e depends on our goals of the computation. 



3.3. 2D system in a disk il. To further verify the effectiveness of 
the penalty method we now test the results in a different domain. We 
now choose a disk Q = {{x, y)\x'^ + 2/^ < !}• The 2D heat equations in 
the polar coordinates x = rcos{6), y = rsin{9) where < ^ < 27r, < 
r < 1 read 



(3.2) 



du Ur uee. , 

— - iy{Urr + — + — ) = /, 

at r 

u\t=o = Mo, 

u\r=i = g. 



Consider the 2D system (13.21) . where < t < 1, and set u = 0.2, / = 
0,(7 = and uo{x,y) = xy, so that g{0) ^ uqIsu- We also set e = 0.1 
the same as before. In this case, we face the singularities almost every- 
where along the unit circle except at the points where x = or y = 0. 
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Section of the solution at t=0 Section of the solution at t=0.02 




Figure 1 1 . The section of the solution aX, 9 = j when 
t is close to 0. 

The effectiveness of this method will be verified with the following nu- 
merical results. 



We first compute the solution of (13. 2p without applying the penalty 
method. The solution is plotted in Fig. [TUl (a) is the graph of the 
solution at t = 0, (b) is the graph of the solution at t = 0.5, and (c) 
is the graph of the solution at t = 1. As we did for the system (13. ip 
for the square, we plot in Fig. [11] the sections (r G (0, 1),^^ = |) of 
the solution at times close to 0. It is clear that, at t = 0, the graph 
displays a sharp gradient around the corner of the time-space axis due 
to the discrepancy between the initial and boundary conditions there, 
and as time evolves, the gradient becomes smoother and smoother. 



To study the error of the system in the disk f2, we define the max- 
imum comparative errors as for the square Q. Hence we plot the L°° 
errors for the 2D system for the disk Q on Fig. [T^l graph (b) is the 
maximum comparative error along the whole time period if we apply 
the penalty method. Because the discrepancy happens at the time- 
space corner, we zoom into the left corner of graph (b) and compare it 
with the error when we do not apply the penalty method. From graph 
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The Error in Both Methods in Short Time Period The Error Along the Whole Time Period 




Figure 12. The maximum comparative error for Ordi- 
nary Finite Element and penalty method in L°° sense 
with e = 0.1 

(a), we observe that the magnitude of the errors at the time-space cor- 
ner are reduced by a factor of 10 if we apply the penalty method. 

In FigHni we plot the maximum comparative error for (13. 2p with a 
fixed mesh at both initial and final steps. At the initial step, the error 
decreases sharply as e increases and remains close to 0, and then it 
becomes flat. At the final step, the error increases almost linearly as 
e increases. The observation also leads to the following conclusion: at 
about e = 0.1, the initial error is minimized while the error at final 
step is well controlled. 

As for the previous example, we shall now look at how the singu- 
larity, induced by the compatibility between the initial and boundary 
data, affects the convergence rates of the numerical scheme. In Fig. 
[13] we plot the maximum errors, at the initial and final time steps, 
with and without the penalty method, against the spatial resolution 
in the loglog scale. We see in Fig. [H] (a) that, without the penalty 
method, the maximum errors do not decrease as the grid refines, which 
demonstrates that the singularity in the solution at the initial time 
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Errors at Initial Step _3 Errors at Final Step 

X 10 




Epsilon Epsilon 
w (b) 

Figure 13. The maximum comparative errors for 
penalty method at both initial and final steps as e vari- 
ants with mesh Ar = ^, A^^ = ^, AT = 

step is serious. With the penalty method, the maximum errors decay 
at roughly the second order. Fig. [T3](b) shows that, with and without 
applying penalty method, the maximum errors at the final steps (t=l) 
decay at approximately the second order. 

3.4. Implementation in a ID System. As we said in the Introduc- 
tion the penalty method applies without any restriction on space di- 
mension. However a number of methods have previously been proposed 
which only apply to space dimension one. Our aim is now to compare 
the efficiency of the penalty method with some of the earlier methods; 
and therefore we can only consider the case of space dimension 1. More 
precisely we will consider the Corrector Methods as proposed in [8]-[T0] 
and compare them with the penalty method for the ID system 



(3.3) 



Ut — vuxx = 0, 0<a;<l, 0<t<l, 
m(x, 0) = Uq 

«(0,t) = (7i(t), u{l,t) = g2{t). 
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Maximum Error at the Initial Time Steps Maximum Error at the Final Time Steps 




-5 5 ' ' 

'1 1.5 2 ■ 1 1.5 2 

logarithm of the number of the logarithm of the number of the 

(a) segment in r (base 1 0) (b) segment in r (base 1 0) 

Figure 14. Decay of the maximum errors. When we 
apply the penalty method here e = 0.1. (a) at the initial 
steps, (b) at the final steps 

Here we set uq{x) = sin{—x + —), giit) = 0, g2(t) = 0, u = 0.2. For 

the Penalty Method, we also set e = 0.1, and for the Corrector Method, 
we have the following choice of correctors offering increasing 

accuracy: 




(3.4) S = { aoSo, (Procedure 1) 

(Procedure 2) 

where ao = 5'i(0)-mo(0), ai = git{0)-Uo^^{0), Sq = / e ^-^ds 



X 



erfc{^^) and Si = I So{x,T)dT. Here Procedure 1 absorbs the O*'* 







order incompatibility {gi{0) ^ uo{0)), and Procedure 2 absorbs both 
the O*'* and 1*** order incompatibilities {gi{0) ^ mo(0) and (7it(0) ^ 
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Without Applying Any Mettiods 




Corrector Method vs Penalty Method 



(a) 




(b) 



Figure 15. Comparative error of the two methods in 
ID system in L°° sense at e = 0.1 



Let u = V + S; we see that v is the solution of the following equation 
Vt — l^Vxx = 0, 0<x<l, 0<t<l, 
(3.5) <^ f (x, 0) = Mo(x), 

viO, t) = g^it) - S{Q, t), v{l, t) = g^it) - S{1, t). 



We choose to solve equation (13. 5p by finite differences. Fig. [15] 
gives the comparison between different methods (Penalty Method and 
Correction Method). Figure [15] (a) gives us the maximum compara- 
tive error of system (13. 3p without applying any methods. Figure [15] 
(b) compares the two methods, zooming into the corner of the time- 
space domain where errors are the largest due to the incompatibility at 
if: = 0. As expected Procedure 2 gives slightly better results than Pro- 
cedure 1. Also the errors with the penalty method are larger than with 
both procedures, but still of comparable magnitude whereas the errors 
without any procedure reach a pick about 6 times larger (4.8 x 10^^ vs 
0.8 X 10^'^). Now we want to vary e in this ID system, FigJT6] shows 
that if e is too small as compared to the mesh, the Penalty Method 
would not reduce the errors at the spatio-temporal corner, but if it is 
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^,Q-3 Errors as epsilon is variant 

'Hf ' ' ' ' ' ' 





Figure 16. The maximum comparative errors for the 
ID system in L°° sense along the time 

an appropriate small number, it could really reduce the errors by more 
than 80%. 

4. Conclusion 

The penalty method gives a way to solve the higher dimensional in- 
compatibility problems. As expected, there exists a solution for system 
(11. 2p which is continuous over [to, 7"], for all to > 0. 

The discrepancy occurs at the time-space corner; we are effectively 
interested in the errors for the initial short time period. The numerical 
simulations for the system with both a square Q and a disk Q yield 
similar results. At the spatio-temporal corner, the magnitudes of the 
errors are reduced by about one order of magnitude by the penalty 
method. Tests are also conducted to study the effects of different val- 
ues of e, the key parameter in the penalty method. We find that with 
an appropriate small value for e, the initial error can be minimized 
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while the error at final step is under well controlled. 

Finally, in space dimension one, when both methods are available 
(penalty method and correction procedures 1 and 2), the penalty method 
gives a slightly larger error than the Correction Procedures 1 and 2; 
but the order of magnitude of the errors are comparable and they are 
all significantly smaller than the errors appearing when no correction 
procedure is implemented . 
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Appendix: The user guide 

The aim is to address the incompatibility issue for the multi-dimensional 
time-dependent linear parabolic equation 

{ut - uAu = /, X eVicR'^, teR+, 
u\t=o = Uo, 
u\dn = g- 

where uqIqq ^ fi'|t=o- So we consider new system instead, namely, for 
£ > fixed. 



{u\ - z/Am^ = /, X G C i?"*, t G 



(43^ ikt + ^{k'-g)=0, teR+, 

( k'{0) = uo\dn- 

We consider for instance the rectangle 0<a;<l,0<y<l and < 
t < 1. We consider the discretization meshes Ax = 1/M, Ay = 1/N 
and At = 1/T, where M, N, T are integers. We use an explicit scheme 
to compute the numerical solution of the original system (14. ip and of 
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the modified system (14.21) . (14.31) . that is respectively: 
(4.4) 

n+l _ n n , n _ n n n , n _ o n 

At ^ Ax2 ^ Ay^ ' 

/or 1 < i < iV - 1, 1 < j < M - 1, 1 < n < T, 

On = ^ij(nAt)|an, /or i = 0,A^ or j = 0, M, 

= uo(iAa;, jAy), /or < z < iV, < j < M. 

for dH]), and , for (g^D-dlJ]): 

At ^ Ax2 A?/2 ^ 

/or 1 < i < iV - 1, 1 < J < A'/ - 1, 1 < n < T, 

>n = kf^U, for ^ = 0, iV or J = 0, M, 
ul^ = uo{iAx,jAy), for 0<i<N, < j < M, 

l.en+1 Ten 



— ^ + -{kr, - g^,,{nAt)) = 0, for z = 0, N or j = 0, M, n>l, 
i\t e 



^ kf. = uo{iAx,jAy), for i = 0,N or j = 0, M. 
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